### Political fragmentation over time
### January 30, 2021


# Data preparation --------------------------------------------------------

### load packages
library(foreign)
library(ggplot2)
library(tidyr)
library(dplyr)
library(maptools)
library(gridExtra)
library(ggmap)
library(data.table)
library(raster)
library(rgdal)
library(sp)
library(grid)
library(gridExtra)
library(Matching)
library(stargazer)
library(lmtest)
library(sandwich)
library(xtable)
library(geosphere)
library(readxl)
library(gganimate)
library(gapminder)
library(gifski)
library(png)
library(transformr)
library(sf)
library(spData)
library(spDataLarge)

### empty environment
rm(list=ls())

### set working directory
setwd("/Users/hanslueders/Dropbox (IPL)/Data copy/Analyses/Fragmentation")




# Historical borders ------------------------------------------------------

### write a function
fragmentation <- function(year){
  # load shapefile
  if(year<1300){
    shp <- readOGR(paste("/Fragmentation Abramson 2017 Shape Files copy/",
                         year,"_poly.shp",sep=""))
  }
  if(year>=1300){
    shp <- readOGR(paste("/Fragmentation Abramson 2017 Shape Files copy/poly_",
                         year,".shp",sep=""))    
  }
  
  # implement name changes (from Abramson)
  shp@data$Name<-ifelse(shp@data$Name =="Balearic islands" ,"Balearic Islands",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name =="Albanians" ,"Albania",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name =="Almavorid_4" ,"Almoravid_4",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name =="Almafi","Amalfi",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name =="Bresancon","Besancon",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name == "Stratsbourg","Strasbourg",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name =="Russians" ,"Russians_1",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name =="Correggio" ,"Corregio" ,as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name =="Wurrtemberg","Wurttemberg",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name =="Ertfurt" ,"Erfurt" ,as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name =="Regensberg" ,"Regensburg"  ,as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name =="Ullster","Ulster",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name =="The  Holy Roman Empire","The Holy Roman Empire",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name =="Fatimad_1","Fatimid_1",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name =="Fatimad_2","Fatimid_2",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name =="Fatimad_3","Fatimid_3",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name =="Fatimad_4","Fatimid_4",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name =="Fatimad_5","Fatimid_5",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name =="Fatimad_6","Fatimid_6",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name =="Kevian Rus","Kievan Rus",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name =="Lithuania","Lithuanians",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name =="Polaban Slavs","Polabian Slavs",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name =="Poblabian Slavs","Polabian Slavs",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name =="Pechengens","Pechenegs",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name =="Kingdom of Burgandy","Kingdom of Burgundy",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name =="Leistner","Leinster",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name =="Noresmen","Norsemen",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name =="South_Slavs","South Slavs",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name == "Strachclyde","Strathclyde",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name == "Umayaa_1_4_1","Ummayad_1_4_1",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name == "Umayaad_1_1","Ummayad_1_1",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name == "Umayaad_1_2","Ummayad_1_2",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name == "Umayaad_1_3","Ummayad_1_3",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name == "Umayaad_1_4","Ummayad_1_4",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name == "Umayaad_1_4_1","Ummayad_1_4_1",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name == "Umayaad_1_4_2","Ummayad_1_4_2",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name == "Umayaad_2","Ummayad_2",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name == "Umayaad_3","Ummayad_3",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name == "Umayaad_3_1","Ummayad_3_1",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name == "Umayaad_3_2","Ummayad_3_2",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name == "Umayaad_4","Ummayad_4",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name == "Umayaad_4_1","Ummayad_4_1",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name == "Umayaad_4_1_1","Ummayad_4_1_1",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name == "Ummayaad_4_1_1","Ummayad_4_1_1",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name == "Umayaad_4_1_2","Ummayad_4_1_2",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name == "Umayaad_4_2","Ummayad_4_2",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name == "Umayaad_5","Ummayad_5",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name == "Umayaard_3_1","Ummayad_3_1",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name == "Umayyad_3_2","Ummayad_3_2",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name == "Geona" , "Genoa",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name == "Almovorads", "Almoravids",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name == "Balearic" | shp@data$Name ==         "Balearic islands"  , "Balearic Islands",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name == "Hammadid", "Hammadids",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name == "Norseman", "Norsemen",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name == "Polostk", "Polotsk",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name =="Munstert", "Munster",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name == "Alhomad" | shp@data$Name ==  "Almohad" , "Almohads",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name == "Bresica"  , "Brescia"  ,as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name ==  "Bulgarians"  ,  "Bulgaria"   ,as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name == "Catile" | shp@data$Name == "Castille", "Castile" ,as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name ==    "Cremoa"   ,   "Cremona"    ,as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name ==    "Dauhphine"   ,    "Dauphine"  ,as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name == "Felte and Belluno"  | shp@data$Name == "Feltre amd Belluno",  "Feltre and Belluno"  ,as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name ==   "Franche Compte" ,    "Franche Comte" ,as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name ==   "Kingdom of the Sword" ,     "Knights of the Sword"  ,as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name ==   "Kingdon of Naples"| shp@data$Name == "The Kingdom of Naples" ,  "Kingdom of Naples"  ,as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name ==   "Pistoa"   ,     "Pistoia"  ,as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name ==  "Portgual"     ,    "Portugal"   ,as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name ==   "The Republic of Novgorod"      ,   "Republic of Novgorod"   ,as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name ==    "Trevisto" ,   "Treviso"   ,as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name == "Entense","Estense",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name =="Freising","Friesberg",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name =="GIlan","Gilan",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name =="Granada","Grenada",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name =="Hapsburgs" ,"Habsburgs",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name == "Il Khanid Mongols", "Il-Khanid Mongols",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name == "Knights Hospitaller", "Knights Hopitaller",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name ==  "Luxembourg", "Luxembourgs",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name ==  "Mongol Empires", "Mongols",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name ==  "Mongol Empire", "Mongols",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name ==  "Vulga Bulgars", "Volga Bulgars"  ,as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name ==  "Wittelsbachs", "Wittlesbachs",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name ==  "Wurtemburg", "Wurttemberg",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name ==  "Hainaut", "Hainault",as.character(shp@data$Name))
  
  shp<-subset(shp,!is.na(shp@data$Name))
  shp@data$X<-as.numeric(as.character(shp@data$X))
  shp@data$Y<-as.numeric(as.character(shp@data$Y))
  shp@data$Area<-as.numeric(as.character(shp@data$Area))
  shp@data$logArea<-log(shp@data$Area)
  shp<-subset(shp,shp@data$Area > 6)
  
  # convert to sf file
  shp2 <- st_as_sf(shp)
  
  # adjust coordinate systems
  st_crs(shp2) <- "+proj=aea +lat_1=43 +lat_2=62 +lat_0=30 +lon_0=10 +x_0=0 +y_0=0 +ellps=intl +units=m +no_defs"
  shp_km <- st_transform(shp2, "+proj=utm +zone=32N +datum=WGS84 +units=km")
  
  # determine centroids
  cent <- centroid(shp)
  cent <- SpatialPoints(cent, proj4string=CRS(as.character("+proj=aea +lat_1=43 +lat_2=62 +lat_0=30 +lon_0=10 +x_0=0 +y_0=0 +ellps=intl +units=m +no_defs")))
  
  # convert and create buffer
  cent <- st_as_sf(cent,coords=c("longitude.y","latitude.y"))
  cent_km <- st_transform(cent, "+proj=utm +zone=32N +datum=WGS84 +units=km")
  cent_buffer50 <- st_buffer(cent_km, 50)
  cent_buffer100 <- st_buffer(cent_km, 100)
  cent_buffer250 <- st_buffer(cent_km, 250)
  cent_buffer500 <- st_buffer(cent_km, 500)
  
  # create output file
  out <- data.frame(matrix(nrow=dim(shp@data[1]), ncol=14))
  colnames(out) <- c("Name", "year", "frag50N", "frag50Meanarea", "frag50Medianarea", 
                     "frag100N", "frag100Meanarea", "frag100Medianarea",
                     "frag250N", "frag250Meanarea", "frag250Medianarea",
                     "frag500N", "frag500Meanarea", "frag500Medianarea")
  out[,1] <- as.character(shp@data$Name)
  out[,2] <- year
  
  
  
  # intersection
  for(j in 1:nrow(out)){
    # define the intersections
    int50 <- st_intersection(cent_buffer50[j,], st_buffer(shp_km, 0))
    int100 <- st_intersection(cent_buffer100[j,], st_buffer(shp_km, 0))
    int250 <- st_intersection(cent_buffer250[j,], st_buffer(shp_km, 0))
    int500 <- st_intersection(cent_buffer500[j,], st_buffer(shp_km, 0))
    
    # compute N, mean area, median area
    out[j,3] <- dim(int50)[1]
    out[j,4] <- mean(int50$Area, na.rm=T)
    out[j,5] <- median(int50$Area, na.rm=T)  
    out[j,6] <- dim(int100)[1]
    out[j,7] <- mean(int100$Area, na.rm=T)
    out[j,8] <- median(int100$Area, na.rm=T)
    out[j,9] <- dim(int250)[1]
    out[j,10] <- mean(int250$Area, na.rm=T)
    out[j,11] <- median(int250$Area, na.rm=T)  
    out[j,12] <- dim(int500)[1]
    out[j,13] <- mean(int500$Area, na.rm=T)
    out[j,14] <- median(int500$Area, na.rm=T)    
  }
  
  # output
  return(out)
}


### Let's try this
# start with year = 1100
out1100 <- fragmentation(year = 1100)
fragmentationDATA <- out1100

# add all other years
for(y in seq(1105,1790, 5)){
  print(y)
  fragmentationDATA <- rbind(fragmentationDATA, fragmentation(year = y))
}

# drop repeat observations
fragmentationDATA2 <- subset(fragmentationDATA, !(Name == "Milan" & year == 1130))
fragmentationDATA2 <- subset(fragmentationDATA2, !(Name == "Milan" & year == 1135 & frag50N == 3))
fragmentationDATA2 <- subset(fragmentationDATA2, !(Name == "Milan" & year == 1170 & frag100N == 7))
fragmentationDATA2 <- subset(fragmentationDATA2, !(Name == "The Holy Roman Empire" & year == 1270 & frag50N == 4))

# save
save(fragmentationDATA, fragmentationDATA2, file="fragmentationDATA.RData")
write.dta(fragmentationDATA2, file="fragmentationDATA.dta")




# 2020 borders ------------------------------------------------------------

### load shapefile
europe <- readOGR("/Users/hanslueders/Dropbox (IPL)/Data copy/GIS data/Europe 2020/ne_110m_admin_0_countries.shp")
europe <- subset(europe, europe@data$CONTINENT == "Europe" & europe@data$SOVEREIGNT != "Russia")
europe.fort <- fortify(europe, region = "SOVEREIGNT")

### determine centroids
cent <- centroid(europe)
cent <- SpatialPoints(cent, proj4string=CRS(as.character("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0")))

### convert and create buffer
cent <- st_as_sf(cent,coords=c("longitude.y","latitude.y"))
cent_km <- st_transform(cent, "+proj=utm +zone=42N +datum=WGS84 +units=km")
cent_buffer50 <- st_buffer(cent_km, 50)
cent_buffer100 <- st_buffer(cent_km, 100)
cent_buffer200 <- st_buffer(cent_km, 200)


### output files
# 50km radius
out50 <- matrix(nrow=length(seq(1100,1790,5)), ncol=39)
out50[,1] <- seq(1100,1790,5)

out50.size.mean <- matrix(nrow=length(seq(1100,1790,5)), ncol=39)
out50.size.mean[,1] <- seq(1100,1790,5)

out50.size.median <- matrix(nrow=length(seq(1100,1790,5)), ncol=39)
out50.size.median[,1] <- seq(1100,1790,5)

# 100km radius
out100 <- matrix(nrow=length(seq(1100,1790,5)), ncol=39)
out100[,1] <- seq(1100,1790,5)

out100.size.mean <- matrix(nrow=length(seq(1100,1790,5)), ncol=39)
out100.size.mean[,1] <- seq(1100,1790,5)

out100.size.median <- matrix(nrow=length(seq(1100,1790,5)), ncol=39)
out100.size.median[,1] <- seq(1100,1790,5)

# 200km radius
out200 <- matrix(nrow=length(seq(1100,1790,5)), ncol=39)
out200[,1] <- seq(1100,1790,5)

out200.size.mean <- matrix(nrow=length(seq(1100,1790,5)), ncol=39)
out200.size.mean[,1] <- seq(1100,1790,5)

out200.size.median <- matrix(nrow=length(seq(1100,1790,5)), ncol=39)
out200.size.median[,1] <- seq(1100,1790,5)


### loop
for(i in 1:length(seq(1100,1790,5))){
  # determine year
  year <- seq(1100,1790,5)[i]
  print(year)
  
  # load shapefile
  if(year<1300){
    shp <- readOGR(paste("/Users/hanslueders/Dropbox (IPL)/Data copy/Shape_Files_For_Anna copy/",
                         year,"_poly.shp",sep=""))
  }
  if(year>=1300){
    shp <- readOGR(paste("/Users/hanslueders/Dropbox (IPL)/Data copy/Shape_Files_For_Anna copy/poly_",
                         year,".shp",sep=""))    
  }
  
  # prepare shapefile
  # convert to sf file
  shp <- st_as_sf(shp)
  
  # adjust coordinate systems
  st_crs(shp) <- "+proj=aea +lat_1=43 +lat_2=62 +lat_0=30 +lon_0=10 +x_0=0 +y_0=0 +ellps=intl +units=m +no_defs"
  shp_km <- st_transform(shp, "+proj=utm +zone=42N +datum=WGS84 +units=km")
  
  
  # intersection
  for(j in 1:dim(cent_buffer50)[1]){
    out50[i,j+1] <- dim(st_intersection(cent_buffer50[j,], st_buffer(shp_km, 0)))[1]
    out100[i,j+1] <- dim(st_intersection(cent_buffer100[j,], st_buffer(shp_km, 0)))[1]
    out200[i,j+1] <- dim(st_intersection(cent_buffer200[j,], st_buffer(shp_km, 0)))[1]
    
    out50.size.mean[i,j+1] <- mean(st_intersection(cent_buffer50[j,], st_buffer(shp_km, 0))$Area, na.rm=T)
    out50.size.median[i,j+1] <- median(st_intersection(cent_buffer50[j,], st_buffer(shp_km, 0))$Area, na.rm=T)
    
    out100.size.mean[i,j+1] <- mean(st_intersection(cent_buffer100[j,], st_buffer(shp_km, 0))$Area, na.rm=T)
    out100.size.median[i,j+1] <- median(st_intersection(cent_buffer100[j,], st_buffer(shp_km, 0))$Area, na.rm=T) 
    
    out200.size.mean[i,j+1] <- mean(st_intersection(cent_buffer200[j,], st_buffer(shp_km, 0))$Area, na.rm=T)
    out200.size.median[i,j+1] <- median(st_intersection(cent_buffer200[j,], st_buffer(shp_km, 0))$Area, na.rm=T)      
  }
}


### turn into data frames
# 50km radius
out50 <- data.frame(out50)
colnames(out50)[1] <- "year"
colnames(out50)[2:39] <- as.character(europe@data$SOVEREIGNT)

out50.size.mean <- data.frame(out50.size.mean)
colnames(out50.size.mean)[1] <- "year"
colnames(out50.size.mean)[2:39] <- as.character(europe@data$SOVEREIGNT)

out50.size.median <- data.frame(out50.size.median)
colnames(out50.size.median)[1] <- "year"
colnames(out50.size.median)[2:39] <- as.character(europe@data$SOVEREIGNT)

# 100km radius
out100 <- data.frame(out100)
colnames(out100)[1] <- "year"
colnames(out100)[2:39] <- as.character(europe@data$SOVEREIGNT)

out100.size.mean <- data.frame(out100.size.mean)
colnames(out100.size.mean)[1] <- "year"
colnames(out100.size.mean)[2:39] <- as.character(europe@data$SOVEREIGNT)

out100.size.median <- data.frame(out100.size.median)
colnames(out100.size.median)[1] <- "year"
colnames(out100.size.median)[2:39] <- as.character(europe@data$SOVEREIGNT)

# 200km radius
out200 <- data.frame(out200)
colnames(out200)[1] <- "year"
colnames(out200)[2:39] <- as.character(europe@data$SOVEREIGNT)

out200.size.mean <- data.frame(out200.size.mean)
colnames(out200.size.mean)[1] <- "year"
colnames(out200.size.mean)[2:39] <- as.character(europe@data$SOVEREIGNT)

out200.size.median <- data.frame(out200.size.median)
colnames(out200.size.median)[1] <- "year"
colnames(out200.size.median)[2:39] <- as.character(europe@data$SOVEREIGNT)


### plot
# 50km radius
p50 <- ggplot() + 
  geom_line(data=out50, aes(x=year, y=Germany, col="Germany", lty="Germany")) +
  geom_line(data=out50, aes(x=year, y=Italy, col="Italy", lty="Italy")) +
  geom_line(data=out50, aes(x=year, y=`United Kingdom`, col="England", lty="England")) +
  geom_line(data=out50, aes(x=year, y=France, col="France", lty="France")) +
  geom_line(data=out50, aes(x=year, y=Spain, col="Spain", lty="Spain")) +  
  scale_x_continuous(breaks=seq(1100,1800,100)) +  
  scale_linetype_manual(values=c(3,4,2,1,5), name="") +
  scale_color_manual(values=c("firebrick3", "steelblue3", "gray25", "forestgreen","gold3"), name="") + 
  xlab("Year") + 
  ylab("Number of states") + 
  theme_classic() + 
  theme(axis.title = element_text(size=12, face="bold"),
        axis.text = element_text(size=10),
        legend.position = "bottom",
        legend.title = element_blank())
p50
ggsave("FragmentationContemporary50km.png", p50, width=8,height=6)

# 100km radius
p100 <- ggplot() + 
  geom_line(data=out100, aes(x=year, y=Germany, col="Germany", lty="Germany")) +
  geom_line(data=out100, aes(x=year, y=Italy, col="Italy", lty="Italy")) +
  geom_line(data=out100, aes(x=year, y=`United Kingdom`, col="England", lty="England")) +
  geom_line(data=out100, aes(x=year, y=France, col="France", lty="France")) +
  geom_line(data=out100, aes(x=year, y=Spain, col="Spain", lty="Spain")) +    
  scale_x_continuous(breaks=seq(1100,1800,100)) +  
  scale_linetype_manual(values=c(3,4,2,1,5), name="") +
  scale_color_manual(values=c("firebrick3", "steelblue3", "gray25", "forestgreen","gold3"), name="") + 
  xlab("Year") + 
  ylab("Number of states") + 
  theme_classic() + 
  theme(axis.title = element_text(size=12, face="bold"),
        axis.text = element_text(size=10),
        legend.position = "bottom",
        legend.title = element_blank())
p100
ggsave("FragmentationContemporary100km.png", p100, width=8,height=6)

# 200km radius
p200 <- ggplot() + 
  geom_line(data=out200, aes(x=year, y=Germany, col="Germany", lty="Germany")) +
  geom_line(data=out200, aes(x=year, y=Italy, col="Italy", lty="Italy")) +
  geom_line(data=out200, aes(x=year, y=`United Kingdom`, col="England", lty="England")) +
  geom_line(data=out200, aes(x=year, y=France, col="France", lty="France")) +
  geom_line(data=out200, aes(x=year, y=Spain, col="Spain", lty="Spain")) +      
  scale_x_continuous(breaks=seq(1100,1800,100)) +
  scale_linetype_manual(values=c(3,4,2,1,5), name="") +
  scale_color_manual(values=c("firebrick3", "steelblue3", "gray25", "forestgreen","gold3"), name="") + 
  xlab("Year") + 
  ylab("Number of states") + 
  theme_classic() + 
  theme(axis.title = element_text(size=12, face="bold"),
        axis.text = element_text(size=10),
        legend.position = "bottom",
        legend.title = element_blank())
p200
ggsave("FragmentationContemporary200km.png", p200, width=8,height=6)


# 200km radius--smoothed
p200 <- ggplot() + 
  geom_smooth(data=out200, aes(x=year, y=Germany, col="Germany", lty="Germany"), se=F, span=.1) +
  geom_smooth(data=out200, aes(x=year, y=Italy, col="Italy", lty="Italy"), se=F, span=.1) +
  geom_smooth(data=out200, aes(x=year, y=`United Kingdom`, col="England", lty="England"), se=F, span=.1) +
  geom_smooth(data=out200, aes(x=year, y=France, col="France", lty="France"), se=F, span=.1) +
  geom_smooth(data=out200, aes(x=year, y=Spain, col="Spain", lty="Spain"), se=F, span=.1) +     
  scale_x_continuous(breaks=seq(1100,1800,100)) +
  scale_linetype_manual(values=c(3,4,2,1,5), name="") +
  scale_color_manual(values=c("firebrick3", "steelblue3", "gray25", "forestgreen","gold3"), name="") + 
  xlab("Year") + 
  ylab("Number of states") + 
  theme_classic() + 
  theme(axis.title = element_text(size=12, face="bold"),
        axis.text = element_text(size=10),
        legend.position = "bottom",
        legend.title = element_blank())
p200
ggsave("FragmentationContemporary200km_smoothed.png", p200, width=8,height=6)


# 200km radius--logged
p200 <- ggplot() + 
  geom_line(data=out200, aes(x=year, y=log(Germany), col="Germany", lty="Germany")) +
  geom_line(data=out200, aes(x=year, y=log(Italy), col="Italy", lty="Italy")) +
  geom_line(data=out200, aes(x=year, y=log(`United Kingdom`), col="England", lty="England")) +
  geom_line(data=out200, aes(x=year, y=log(France), col="France", lty="France")) +
  geom_line(data=out200, aes(x=year, y=log(Spain), col="Spain", lty="Spain")) +    
  scale_x_continuous(breaks=seq(1100,1800,100)) +
  scale_linetype_manual(values=c(3,4,2,1,5), name="") +
  scale_color_manual(values=c("firebrick3", "steelblue3", "gray25", "forestgreen","gold3"), name="") + 
  xlab("Year") + 
  ylab("Number of states (log)") + 
  theme_classic() + 
  theme(axis.title = element_text(size=12, face="bold"),
        axis.text = element_text(size=10),
        legend.position = "bottom",
        legend.title = element_blank())
p200
ggsave("FragmentationContemporary200km_logged.png", p200, width=8,height=6)

# 200km radius--mean area (log)
p200 <- ggplot() + 
  geom_line(data=out200.size.mean, aes(x=year, y=log(Germany), col="Germany", lty="Germany")) +
  geom_line(data=out200.size.mean, aes(x=year, y=log(Italy), col="Italy", lty="Italy")) +
  geom_line(data=out200.size.mean, aes(x=year, y=log(`United Kingdom`), col="England", lty="England")) +
  geom_line(data=out200.size.mean, aes(x=year, y=log(France), col="France", lty="France")) +
  geom_line(data=out200.size.mean, aes(x=year, y=log(Spain), col="Spain", lty="Spain")) +  
  scale_x_continuous(breaks=seq(1100,1800,100)) +
  scale_linetype_manual(values=c(3,4,2,1,5), name="") +
  scale_color_manual(values=c("firebrick3", "steelblue3", "gray25", "forestgreen","gold3"), name="") + 
  xlab("Year") + 
  ylab("Mean area (log km2)") + 
  theme_classic() + 
  theme(axis.title = element_text(size=12, face="bold"),
        axis.text = element_text(size=10),
        legend.position = "bottom",
        legend.title = element_blank())
p200
ggsave("FragmentationContemporary200km_meanarea.png", p200, width=8,height=6)


# 200km radius--median area (log)
p200 <- ggplot() + 
  geom_line(data=out200.size.median, aes(x=year, y=log(Germany), col="Germany", lty="Germany")) +
  geom_line(data=out200.size.median, aes(x=year, y=log(Italy), col="Italy", lty="Italy")) +
  geom_line(data=out200.size.median, aes(x=year, y=log(`United Kingdom`), col="England", lty="England")) +
  geom_line(data=out200.size.median, aes(x=year, y=log(France), col="France", lty="France")) +
  geom_line(data=out200.size.median, aes(x=year, y=log(Spain), col="Spain", lty="Spain")) +  
  scale_x_continuous(breaks=seq(1100,1800,100)) +
  scale_linetype_manual(values=c(3,4,2,1,5), name="") +
  scale_color_manual(values=c("firebrick3", "steelblue3", "gray25", "forestgreen","gold3"), name="") + 
  xlab("Year") + 
  ylab("Median area (log km2)") + 
  theme_classic() + 
  theme(axis.title = element_text(size=12, face="bold"),
        axis.text = element_text(size=10),
        legend.position = "bottom",
        legend.title = element_blank())
p200
ggsave("FragmentationContemporary200km_medianarea.png", p200, width=8,height=6)





# Number of states within 1900 boundaries ---------------------------------

### load shapefile
A <- readOGR("/Users/hanslueders/Dropbox (IPL)/Data copy/GIS data/Europe Historical/01 Europe Main/Europe_1900_v.1.1.shp")
#A.fort <- fortify(A, region="COUNTRY")

### convert to sf
A2 <- st_as_sf(A,coords=c("longitude.y","latitude.y"))
st_crs(A2) <- "+proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs"

### output
out <- matrix(nrow=length(seq(1100,1790,5)), ncol=7)
out[,1] <- seq(1100,1790,5)

out.size.mean <- matrix(nrow=length(seq(1100,1790,5)), ncol=7)
out.size.mean[,1] <- seq(1100,1790,5)

out.size.median <- matrix(nrow=length(seq(1100,1790,5)), ncol=7)
out.size.median[,1] <- seq(1100,1790,5)

### loop
for(i in 1:length(seq(1100,1790,5))){
  # determine year
  year <- seq(1100,1790,5)[i]
  print(year)
  
  # load shapefile
  if(year<1300){
    shp <- readOGR(paste("/Users/hanslueders/Dropbox (IPL)/Data copy/Shape_Files_For_Anna copy/",
                         year,"_poly.shp",sep=""))
  }
  if(year>=1300){
    shp <- readOGR(paste("/Users/hanslueders/Dropbox (IPL)/Data copy/Shape_Files_For_Anna copy/poly_",
                         year,".shp",sep=""))    
  }
  
  # adjust CRS
  shp2 <- st_as_sf(shp)
  
  # adjust CRS
  st_crs(shp2) <- "+proj=aea +lat_1=43 +lat_2=62 +lat_0=30 +lon_0=10 +x_0=0 +y_0=0 +ellps=intl +units=m +no_defs"
  shp2 <- st_transform(shp2, "+proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs")
  
  
  # determine overlap
  out[i,2] <- length(unique(unlist(st_intersects(A2[A2$COUNTRY==50,],shp2)))) # England
  out[i,3] <- length(unique(unlist(st_intersects(A2[A2$COUNTRY==80,],shp2)))) # Germany
  out[i,4] <- length(unique(unlist(st_intersects(A2[A2$COUNTRY==110,],shp2)))) # Italy
  out[i,5] <- length(unique(unlist(st_intersects(A2[A2$COUNTRY != 80 & A2$COUNTRY != 110,],shp2)))) # neither Germany nor Italy
  out[i,6] <- length(unique(unlist(st_intersects(A2[A2$COUNTRY==181,],shp2)))) # Spain
  out[i,7] <- length(unique(unlist(st_intersects(A2[A2$COUNTRY==70,],shp2)))) # France
  
  # determine size
  out.size.mean[i,2] <- mean(shp2$Area[unique(unlist(st_intersects(A2[A2$COUNTRY==50,],shp2)))], na.rm=T) # England
  out.size.mean[i,3] <- mean(shp2$Area[unique(unlist(st_intersects(A2[A2$COUNTRY==80,],shp2)))], na.rm=T) # Germany
  out.size.mean[i,4] <- mean(shp2$Area[unique(unlist(st_intersects(A2[A2$COUNTRY==110,],shp2)))], na.rm=T) # Italy
  out.size.mean[i,5] <- mean(shp2$Area[unique(unlist(st_intersects(A2[A2$COUNTRY != 80 & A2$COUNTRY != 110,],shp2)))], na.rm=T) # neither Germany nor Italy  
  out.size.mean[i,6] <- mean(shp2$Area[unique(unlist(st_intersects(A2[A2$COUNTRY==181,],shp2)))], na.rm=T) # Spain
  out.size.mean[i,7] <- mean(shp2$Area[unique(unlist(st_intersects(A2[A2$COUNTRY==70,],shp2)))], na.rm=T) # France
  
  out.size.median[i,2] <- median(shp2$Area[unique(unlist(st_intersects(A2[A2$COUNTRY==50,],shp2)))], na.rm=T) # England
  out.size.median[i,3] <- median(shp2$Area[unique(unlist(st_intersects(A2[A2$COUNTRY==80,],shp2)))], na.rm=T) # Germany
  out.size.median[i,4] <- median(shp2$Area[unique(unlist(st_intersects(A2[A2$COUNTRY==110,],shp2)))], na.rm=T) # Italy
  out.size.median[i,5] <- median(shp2$Area[unique(unlist(st_intersects(A2[A2$COUNTRY != 80 & A2$COUNTRY != 110,],shp2)))], na.rm=T) # neither Germany nor Italy  
  out.size.median[i,6] <- median(shp2$Area[unique(unlist(st_intersects(A2[A2$COUNTRY==181,],shp2)))], na.rm=T) # Spain
  out.size.median[i,7] <- median(shp2$Area[unique(unlist(st_intersects(A2[A2$COUNTRY==70,],shp2)))], na.rm=T) # France
  
}


### Manually add the data for the 19th century
# Germany
germany <- cbind(c(1808,1848,1867,1871,1900), rep(1,5), c(40,40,22,1,1), rep(NA,5), rep(NA,5),rep(NA,5),rep(NA,5))

# Italy
italy <- cbind(c(1828,1829,1847,1857,1860,1861,1866,1870,1900), rep(1,9), rep(NA, 9), c(11,10,9,7,4,4,2,1,1), rep(NA, 9), rep(NA, 9), rep(NA, 9))

# France
france <- cbind(c(1850,1870,1900), rep(1,3), rep(NA,3), rep(NA,3), rep(NA,3), rep(NA,3), c(2,1,1))

# Spain
spain <- cbind(c(1790,1850,1900), rep(1,3), rep(NA,3), rep(NA,3), rep(NA,3), c(3,1,1), rep(NA,3))


# merge to master file
out2 <- rbind(out, germany)
out2 <- rbind(out2, italy)
out2 <- rbind(out2, france)
out2 <- rbind(out2, spain)



### convert to data frame
out <- data.frame(out)
colnames(out) <- c("year", "England","Germany","Italy","RestEurope", "Spain", "France")

out2 <- data.frame(out2)
colnames(out2) <- c("year", "England","Germany","Italy","RestEurope", "Spain", "France")

out.size.mean <- data.frame(out.size.mean)
colnames(out.size.mean) <- c("year", "England","Germany","Italy","RestEurope", "Spain", "France")

out.size.median <- data.frame(out.size.median)
colnames(out.size.median) <- c("year", "England","Germany","Italy","RestEurope", "Spain", "France")




### save the data
save(out, out2, out.size.mean, out.size.median,
     file = "FragmentationData1900borders.RData")
#load("FragmentationData1900borders.RData")


### Plots
# number of states
p <- ggplot() + 
  geom_line(data=out, aes(x=year, y=England, col="England", lty="England"),size=0.8) +
  geom_line(data=out, aes(x=year, y=Germany, col="Germany", lty="Germany"),size=0.8) +
  geom_line(data=out, aes(x=year, y=Italy, col="Italy", lty="Italy"),size=0.8) +
  geom_line(data=out, aes(x=year, y=France, col="France", lty="France"),size=0.8) +
  geom_line(data=out, aes(x=year, y=Spain, col="Spain", lty="Spain"),size=0.8) +
  scale_x_continuous(breaks=seq(1100,1800,100)) +
  scale_linetype_manual(values=c(3,4,2,1,5), name="") +
  scale_color_manual(values=c("gray25", "gray25", "gray50", "gray50","gray75"), name="") +  
  #scale_color_manual(values=c("firebrick3", "steelblue3", "gray25", "forestgreen","gold3"), name="") +  
  xlab("Year") + 
  ylab("Number of states") + 
  theme_classic() + 
  theme(axis.title.x = element_blank(),
        axis.title.y = element_text(size=16, face="bold"),
        axis.text.x = element_text(size=14),
        axis.text.y = element_text(size=14),
        legend.position = "bottom",
        legend.title = element_blank(),
        legend.text = element_text(size=12))
p
ggsave("Fragmentation1900borders.png", p, width=8,height=5)

# number of states NEW
p <- ggplot() + 
  geom_line(data=out2[!is.na(out2$England),], aes(x=year, y=England, col="England", lty="England"), size=0.8) +
  geom_line(data=out2[!is.na(out2$Germany),], aes(x=year, y=Germany, col="Germany", lty="Germany"), size=0.8) +
  geom_line(data=out2[!is.na(out2$Italy),], aes(x=year, y=Italy, col="Italy", lty="Italy"), size=0.8) +
  geom_line(data=out2[!is.na(out2$France),], aes(x=year, y=France, col="France", lty="France"), size=0.8) +
  geom_line(data=out2[!is.na(out2$Spain),], aes(x=year, y=Spain, col="Spain", lty="Spain"), size=0.8) +
  scale_x_continuous(breaks=seq(1100,1900,100)) +
  scale_linetype_manual(values=c(3,4,2,1,5), name="") +
  scale_color_manual(values=c("gray25", "gray25", "gray50", "gray50","gray75"), name="") +  
  #scale_color_manual(values=c("firebrick3", "steelblue3", "gray25", "forestgreen","gold3"), name="") +  
  xlab("Year") + 
  ylab("Number of states") + 
  theme_classic() + 
  theme(axis.title.x = element_blank(),
        axis.title.y = element_text(size=16, face="bold"),
        axis.text.x = element_text(size=14),
        axis.text.y = element_text(size=14),
        legend.position = "bottom",
        legend.title = element_blank(),
        legend.text = element_text(size=12))
p
ggsave("Fragmentation1900borders_v2.png", p, width=8,height=5)
ggsave("Fragmentation1900borders_v2.pdf", p, width=8,height=5)

# number of states NEW in color
p <- ggplot() + 
  geom_line(data=out2[!is.na(out2$England) & out2$year < 1800,], aes(x=year, y=England, col="England", lty="England"), size=0.8) +
  geom_line(data=out2[!is.na(out2$Germany) & out2$year < 1800,], aes(x=year, y=Germany, col="Germany", lty="Germany"), size=0.8) +
  geom_line(data=out2[!is.na(out2$Italy) & out2$year < 1800,], aes(x=year, y=Italy, col="Italy", lty="Italy"), size=0.8) +
  geom_line(data=out2[!is.na(out2$France) & out2$year < 1800,], aes(x=year, y=France, col="France", lty="France"), size=0.8) +
  geom_line(data=out2[!is.na(out2$Spain) & out2$year < 1800,], aes(x=year, y=Spain, col="Spain", lty="Spain"), size=0.8) +
  scale_x_continuous(breaks=seq(1100,1800,100)) +
  scale_linetype_manual(values=c(3,4,2,1,5), name="") +
  #scale_color_manual(values=c("gray25", "gray25", "gray50", "gray50","gray75"), name="") +  
  scale_color_manual(values=c("firebrick3", "steelblue3", "gray25", "forestgreen","gold3"), name="") +  
  xlab("Year") + 
  ylab("Number of states") + 
  theme_classic() + 
  theme(axis.title.x = element_blank(),
        axis.title.y = element_text(size=16, face="bold"),
        axis.text.x = element_text(size=14),
        axis.text.y = element_text(size=14),
        legend.position = "bottom",
        legend.title = element_blank(),
        legend.text = element_text(size=12))
p
ggsave("Fragmentation1900borders_v2col.png", p, width=8,height=5)
ggsave("Fragmentation1900borders_v2col.pdf", p, width=8,height=5)

# number of states (log)
p <- ggplot() + 
  geom_line(data=out, aes(x=year, y=log(England), col="England", lty="England"), size=0.8) +
  geom_line(data=out, aes(x=year, y=log(Germany), col="Germany", lty="Germany"), size=0.8) +
  geom_line(data=out, aes(x=year, y=log(Italy), col="Italy", lty="Italy"), size=0.8) +
  geom_line(data=out, aes(x=year, y=log(France), col="France", lty="France"), size=0.8) +
  geom_line(data=out, aes(x=year, y=log(Spain), col="Spain", lty="Spain"), size=0.8) +
  scale_x_continuous(breaks=seq(1100,1800,100)) +
  scale_y_continuous(breaks=seq(0,5,1)) +  
  scale_linetype_manual(values=c(3,4,2,1,5), name="") +
  scale_color_manual(values=c("gray25", "gray25", "gray50", "gray50","gray75"), name="") +  
  #scale_color_manual(values=c("firebrick3", "steelblue3", "gray25", "forestgreen","gold3"), name="") +  
  xlab("Year") + 
  ylab("Number of states (log)") + 
  theme_classic() + 
  theme(axis.title.x = element_blank(),
        axis.title.y = element_text(size=16, face="bold"),
        axis.text.x = element_text(size=14),
        axis.text.y = element_text(size=14),
        legend.position = "bottom",
        legend.title = element_blank(),
        legend.text = element_text(size=12))
p
ggsave("Fragmentation1900borders_logged.png", p, width=8,height=5)

# number of states (log) NEW
p <- ggplot() + 
  geom_line(data=out2[!is.na(out2$England),], aes(x=year, y=log(England), col="England", lty="England"), size=0.8) +
  geom_line(data=out2[!is.na(out2$Germany),], aes(x=year, y=log(Germany), col="Germany", lty="Germany"), size=0.8) +
  geom_line(data=out2[!is.na(out2$Italy),], aes(x=year, y=log(Italy), col="Italy", lty="Italy"), size=0.8) +
  geom_line(data=out2[!is.na(out2$France),], aes(x=year, y=log(France), col="France", lty="France"), size=0.8) +
  geom_line(data=out2[!is.na(out2$Spain),], aes(x=year, y=log(Spain), col="Spain", lty="Spain"), size=0.8) +
  scale_x_continuous(breaks=seq(1100,1900,100)) +
  scale_y_continuous(breaks=seq(0,5,1)) +
  scale_linetype_manual(values=c(3,4,2,1,5), name="") +
  scale_color_manual(values=c("gray25", "gray25", "gray50", "gray50","gray75"), name="") +  
  #scale_color_manual(values=c("firebrick3", "steelblue3", "gray25", "forestgreen","gold3"), name="") +  
  xlab("Year") + 
  ylab("Number of states (log)") + 
  theme_classic() + 
  theme(axis.title.x = element_blank(),
        axis.title.y = element_text(size=16, face="bold"),
        axis.text.x = element_text(size=14),
        axis.text.y = element_text(size=14),
        legend.position = "bottom",
        legend.title = element_blank(),
        legend.text = element_text(size=12))
p
ggsave("Fragmentation1900borders_logged_v2.png", p, width=8,height=5)

# number of states (log) NEW in color
p <- ggplot() + 
  geom_line(data=out2[!is.na(out2$England),], aes(x=year, y=log(England), col="England", lty="England"), size=0.8) +
  geom_line(data=out2[!is.na(out2$Germany),], aes(x=year, y=log(Germany), col="Germany", lty="Germany"), size=0.8) +
  geom_line(data=out2[!is.na(out2$Italy),], aes(x=year, y=log(Italy), col="Italy", lty="Italy"), size=0.8) +
  geom_line(data=out2[!is.na(out2$France),], aes(x=year, y=log(France), col="France", lty="France"), size=0.8) +
  geom_line(data=out2[!is.na(out2$Spain),], aes(x=year, y=log(Spain), col="Spain", lty="Spain"), size=0.8) +
  scale_x_continuous(breaks=seq(1100,1900,100)) +
  scale_y_continuous(breaks=seq(0,5,1)) +
  scale_linetype_manual(values=c(3,4,2,1,5), name="") +
  #scale_color_manual(values=c("gray25", "gray25", "gray50", "gray50","gray75"), name="") +  
  scale_color_manual(values=c("firebrick3", "steelblue3", "gray25", "forestgreen","gold3"), name="") +  
  xlab("Year") + 
  ylab("Number of states (log)") + 
  theme_classic() + 
  theme(axis.title.x = element_blank(),
        axis.title.y = element_text(size=16, face="bold"),
        axis.text.x = element_text(size=14),
        axis.text.y = element_text(size=14),
        legend.position = "bottom",
        legend.title = element_blank(),
        legend.text = element_text(size=12))
p
ggsave("Fragmentation1900borders_logged_v2col.png", p, width=8,height=5)

# number of states (smoothed)
p <- ggplot() + 
  geom_smooth(data=out, aes(x=year, y=England, col="England", lty="England"), se=F, span=.1) +
  geom_smooth(data=out, aes(x=year, y=Germany, col="Germany", lty="Germany"), se=F, span=.1) +
  geom_smooth(data=out, aes(x=year, y=Italy, col="Italy", lty="Italy"), se=F, span=.1) +
  geom_smooth(data=out, aes(x=year, y=France, col="France", lty="France"), se=F, span=.1) +
  geom_smooth(data=out, aes(x=year, y=Spain, col="Spain", lty="Spain"), se=F, span=.1) +  
  scale_x_continuous(breaks=seq(1100,1800,100)) +
  scale_linetype_manual(values=c(3,4,2,1,5), name="") +
  scale_color_manual(values=c("firebrick3", "steelblue3", "gray25", "forestgreen","gold3"), name="") +  
  xlab("Year") + 
  ylab("Number of states") + 
  theme_classic() + 
  theme(axis.title = element_text(size=12, face="bold"),
        axis.text = element_text(size=10),
        legend.position = "bottom",
        legend.title = element_blank())
p
ggsave("Fragmentation1900borders_smoothed.png", p, width=8,height=5)


# number of states (logged mean area)
p <- ggplot() + 
  geom_line(data=out.size.mean, aes(x=year, y=log(England), col="England", lty="England")) +
  geom_line(data=out.size.mean, aes(x=year, y=log(Germany), col="Germany", lty="Germany")) +
  geom_line(data=out.size.mean, aes(x=year, y=log(Italy), col="Italy", lty="Italy")) +
  geom_line(data=out.size.mean, aes(x=year, y=log(France), col="France", lty="France")) +
  geom_line(data=out.size.mean, aes(x=year, y=log(Spain), col="Spain", lty="Spain")) +
  scale_x_continuous(breaks=seq(1100,1800,100)) +
  scale_linetype_manual(values=c(3,4,2,1,5), name="") +
  scale_color_manual(values=c("firebrick3", "steelblue3", "gray25", "forestgreen","gold3"), name="") +    
  xlab("Year") + 
  ylab("Mean area (log km2)") + 
  theme_classic() + 
  theme(axis.title = element_text(size=12, face="bold"),
        axis.text = element_text(size=10),
        legend.position = "bottom",
        legend.title = element_blank())
p
ggsave("Fragmentation1900borders_meanarea.png", p, width=8,height=5)


# number of states (logged median area)
p <- ggplot() + 
  geom_line(data=out.size.median, aes(x=year, y=log(England), col="England", lty="England")) +
  geom_line(data=out.size.median, aes(x=year, y=log(Germany), col="Germany", lty="Germany")) +
  geom_line(data=out.size.median, aes(x=year, y=log(Italy), col="Italy", lty="Italy")) +
  geom_line(data=out.size.median, aes(x=year, y=log(France), col="France", lty="France")) +
  geom_line(data=out.size.median, aes(x=year, y=log(Spain), col="Spain", lty="Spain")) +
  scale_x_continuous(breaks=seq(1100,1800,100)) +
  scale_linetype_manual(values=c(3,4,2,1,5), name="") +
  scale_color_manual(values=c("firebrick3", "steelblue3", "gray25", "forestgreen","gold3"), name="") +    
  xlab("Year") + 
  ylab("Median area (log km2)") + 
  theme_classic() + 
  theme(axis.title = element_text(size=12, face="bold"),
        axis.text = element_text(size=10),
        legend.position = "bottom",
        legend.title = element_blank())
p
ggsave("Fragmentation1900borders_meanarea.png", p, width=8,height=5)






# Fragmentation using 1900 borders ----------------------------------------

### load 1900 shapefile
A <- readOGR("/Users/hanslueders/Dropbox (IPL)/Data copy/GIS data/Europe Historical/01 Europe Main/Europe_1900_v.1.1.shp")

### write a function
fragmentation1900 <- function(year){
  # load shapefile
  if(year<1300){
    shp <- readOGR(paste("/Users/hanslueders/Dropbox (IPL)/Data copy/Shape_Files_For_Anna copy/",
                         year,"_poly.shp",sep=""))
  }
  if(year>=1300){
    shp <- readOGR(paste("/Users/hanslueders/Dropbox (IPL)/Data copy/Shape_Files_For_Anna copy/poly_",
                         year,".shp",sep=""))    
  }
  
  # implement name changes (from Abramson)
  shp@data$Name<-ifelse(shp@data$Name =="Balearic islands" ,"Balearic Islands",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name =="Albanians" ,"Albania",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name =="Almavorid_4" ,"Almoravid_4",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name =="Almafi","Amalfi",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name =="Bresancon","Besancon",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name == "Stratsbourg","Strasbourg",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name =="Russians" ,"Russians_1",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name =="Correggio" ,"Corregio" ,as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name =="Wurrtemberg","Wurttemberg",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name =="Ertfurt" ,"Erfurt" ,as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name =="Regensberg" ,"Regensburg"  ,as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name =="Ullster","Ulster",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name =="The  Holy Roman Empire","The Holy Roman Empire",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name =="Fatimad_1","Fatimid_1",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name =="Fatimad_2","Fatimid_2",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name =="Fatimad_3","Fatimid_3",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name =="Fatimad_4","Fatimid_4",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name =="Fatimad_5","Fatimid_5",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name =="Fatimad_6","Fatimid_6",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name =="Kevian Rus","Kievan Rus",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name =="Lithuania","Lithuanians",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name =="Polaban Slavs","Polabian Slavs",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name =="Poblabian Slavs","Polabian Slavs",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name =="Pechengens","Pechenegs",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name =="Kingdom of Burgandy","Kingdom of Burgundy",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name =="Leistner","Leinster",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name =="Noresmen","Norsemen",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name =="South_Slavs","South Slavs",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name == "Strachclyde","Strathclyde",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name == "Umayaa_1_4_1","Ummayad_1_4_1",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name == "Umayaad_1_1","Ummayad_1_1",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name == "Umayaad_1_2","Ummayad_1_2",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name == "Umayaad_1_3","Ummayad_1_3",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name == "Umayaad_1_4","Ummayad_1_4",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name == "Umayaad_1_4_1","Ummayad_1_4_1",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name == "Umayaad_1_4_2","Ummayad_1_4_2",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name == "Umayaad_2","Ummayad_2",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name == "Umayaad_3","Ummayad_3",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name == "Umayaad_3_1","Ummayad_3_1",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name == "Umayaad_3_2","Ummayad_3_2",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name == "Umayaad_4","Ummayad_4",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name == "Umayaad_4_1","Ummayad_4_1",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name == "Umayaad_4_1_1","Ummayad_4_1_1",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name == "Ummayaad_4_1_1","Ummayad_4_1_1",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name == "Umayaad_4_1_2","Ummayad_4_1_2",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name == "Umayaad_4_2","Ummayad_4_2",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name == "Umayaad_5","Ummayad_5",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name == "Umayaard_3_1","Ummayad_3_1",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name == "Umayyad_3_2","Ummayad_3_2",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name == "Geona" , "Genoa",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name == "Almovorads", "Almoravids",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name == "Balearic" | shp@data$Name ==         "Balearic islands"  , "Balearic Islands",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name == "Hammadid", "Hammadids",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name == "Norseman", "Norsemen",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name == "Polostk", "Polotsk",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name =="Munstert", "Munster",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name == "Alhomad" | shp@data$Name ==  "Almohad" , "Almohads",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name == "Bresica"  , "Brescia"  ,as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name ==  "Bulgarians"  ,  "Bulgaria"   ,as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name == "Catile" | shp@data$Name == "Castille", "Castile" ,as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name ==    "Cremoa"   ,   "Cremona"    ,as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name ==    "Dauhphine"   ,    "Dauphine"  ,as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name == "Felte and Belluno"  | shp@data$Name == "Feltre amd Belluno",  "Feltre and Belluno"  ,as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name ==   "Franche Compte" ,    "Franche Comte" ,as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name ==   "Kingdom of the Sword" ,     "Knights of the Sword"  ,as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name ==   "Kingdon of Naples"| shp@data$Name == "The Kingdom of Naples" ,  "Kingdom of Naples"  ,as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name ==   "Pistoa"   ,     "Pistoia"  ,as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name ==  "Portgual"     ,    "Portugal"   ,as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name ==   "The Republic of Novgorod"      ,   "Republic of Novgorod"   ,as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name ==    "Trevisto" ,   "Treviso"   ,as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name == "Entense","Estense",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name =="Freising","Friesberg",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name =="GIlan","Gilan",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name =="Granada","Grenada",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name =="Hapsburgs" ,"Habsburgs",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name == "Il Khanid Mongols", "Il-Khanid Mongols",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name == "Knights Hospitaller", "Knights Hopitaller",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name ==  "Luxembourg", "Luxembourgs",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name ==  "Mongol Empires", "Mongols",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name ==  "Mongol Empire", "Mongols",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name ==  "Vulga Bulgars", "Volga Bulgars"  ,as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name ==  "Wittelsbachs", "Wittlesbachs",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name ==  "Wurtemburg", "Wurttemberg",as.character(shp@data$Name))
  shp@data$Name<-ifelse(shp@data$Name ==  "Hainaut", "Hainault",as.character(shp@data$Name))
  
  shp<-subset(shp,!is.na(shp@data$Name))
  shp@data$X<-as.numeric(as.character(shp@data$X))
  shp@data$Y<-as.numeric(as.character(shp@data$Y))
  shp@data$Area<-as.numeric(as.character(shp@data$Area))
  shp@data$logArea<-log(shp@data$Area)
  shp<-subset(shp,shp@data$Area > 6)
  
  # determine centroids
  cent <- centroid(shp)
  cent <- SpatialPoints(cent, proj4string=CRS(as.character("+proj=aea +lat_1=43 +lat_2=62 +lat_0=30 +lon_0=10 +x_0=0 +y_0=0 +ellps=intl +units=m +no_defs")))
  
  # adjust CRS
  cent2 <- spTransform(cent, CRS("+proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs"))
  
  # assign to 1900 state borders
  by1900 <- over(cent2, A)  
  shp@data$country1900 <- by1900$COUNTRY
  
  # collapse by state
  by1900$count <- 1
  by1900.agg <- aggregate(by1900$count, by=list(country1900 = by1900$COUNTRY), FUN=sum, na.rm=T)
  colnames(by1900.agg) <- c("country1900", "fragN1900")
  
  # merge
  shp@data <- merge(shp@data, by1900.agg, by="country1900", all.x=T)
  
  # select data
  out <- subset(shp@data, select = c(Name, fragN1900))
  out$year <- year

  # output
  return(out)
}


### Implement the function
# start with year = 1100
out1100 <- fragmentation1900(year = 1100)
fragmentationDATA1900 <- out1100

# add all other years
for(y in seq(1105,1790, 5)){
  print(y)
  fragmentationDATA1900 <- rbind(fragmentationDATA1900, fragmentation1900(year = y))
}


### drop repeat observations
fragmentationDATA1900$fragN1900[fragmentationDATA1900$Name == "Milan" & fragmentationDATA1900$year == 1130] <- 8
fragmentationDATA1900$fragN1900[fragmentationDATA1900$Name == "Milan" & fragmentationDATA1900$year == 1135] <- 8
fragmentationDATA1900$fragN1900[fragmentationDATA1900$Name == "The Holy Roman Empire" & fragmentationDATA1900$year == 1270] <- 135
fragmentationDATA1900 <- fragmentationDATA1900[duplicated(fragmentationDATA1900) == F,]

fragmentationDATA2 <- subset(fragmentationDATA2, !(Name == "Milan" & year == 1135 & frag50N == 3))
fragmentationDATA2 <- subset(fragmentationDATA2, !(Name == "Milan" & year == 1170 & frag100N == 7))
fragmentationDATA2 <- subset(fragmentationDATA2, !(Name == "The Holy Roman Empire" & year == 1270 & frag50N == 4))


### save
save(fragmentationDATA1900, file="fragmentationDATA1900.RData")
write.dta(fragmentationDATA1900, file="fragmentationDATA1900.dta")


